Warming waters lead to increased habitat suitability for juvenile bull sharks (Carcharhinus leucas)

Coastal ecosystems are highly vulnerable to the impacts of climate change and other stressors, including urbanization and overfishing. Consequently, distributions of coastal fish have begun to change, particularly in response to increasing temperatures linked to climate change. However, few studies have evaluated how natural and anthropogenic disturbances can alter species distributions in conjunction with geophysical habitat alterations, such as changes to land use and land cover (LU/LC). Here, we examine the spatiotemporal changes in the distribution of juvenile bull sharks (Carcharhinus leucas) using a multi-decadal fishery-independent survey of coastal Alabama. Using a boosted regression tree (BRT) modeling framework, we assess the covariance of environmental conditions (sea surface temperature, depth, salinity, dissolved oxygen, riverine discharge, Chl-a) as well as historic changes to LU/LC to the distribution of bull sharks. Species distribution models resultant from BRTs for early (2003–2005) and recent (2018–2020) monitoring periods indicated a mean increase in habitat suitability (i.e., probability of capture) for juvenile bull sharks from 0.028 to 0.082, concomitant with substantial increases in mean annual temperature (0.058°C/yr), Chl-a (2.32 mg/m3), and urbanization (increased LU/LC) since 2000. These results align with observed five-fold increases in the relative abundance of juvenile bull sharks across the study period and demonstrate the impacts of changing environmental conditions on their distribution and relative abundance. As climate change persists, coastal communities will continue to change, altering the structure of ecological communities and the success of nearshore fisheries.


Study site
Mobile Bay is a shallow, drowned river valley estuary that is highly variable along the transitional river-sea gradient.It is characterized by an extensive delta region at its northern extent, where the Mobile River empties 90% of the freshwater input into the system 36 .Strong seasonal inflow drives salinity trends, with high flow, low salinity periods occurring February-April and low flow, high salinity occurring August-October 36 .River discharge entering the northern bay has a range exceeding two orders of magnitude, making it an extension of the river during peak flows in winter and spring 37 .The bay is moderately stratified year-round, with the strongest gradients occurring in the springtime as a response to increased river discharge 36 .An overall latitudinal salinity gradient is apparent, with fresher waters at the northern portion of Mobile Bay becoming increasingly saline and stable southward, where Gulf of Mexico waters are exchanged at Main Pass 36 .The southern extent of the bay exhibits less thermal range, where it tends to be warmer in the winter and colder in the summer than the northern portion of the upper bay 36 .The Coriolis force and the triangular shoreline configuration deflect Gulf waters along the eastern shore of Mobile Bay, while fresh waters are deflected to the western shore 36 .Overall, thermal stratification is minimal, with less than 1 °C difference between surface and bottom level waters 36 .

Field work
Year-round gillnet surveys have been conducted by the Alabama Department of Conservation and Natural Resources-Marine Resources Division (AMRD) since 2001.For this survey, small and large mesh gillnets were used.Small mesh gillnets were divided into five panels (2.44 × 45.7 m), with mesh sizes ranging from 5.1 to 10.2 cm, each panel increasing incrementally by 1.27 cm.Environmental conditions permitting, these nets were set perpendicular to the shoreline with the smallest mesh size closest to shore.Large mesh gillnets were divided into four panels (2.44 × 45.7 m) ranging from 11.4 to 15.2 cm, with each panel increasing incrementally by 1.27 cm.Large mesh nets were set parallel to the shore to a maximum depth of 8 ft.Both net types were soaked for one hour, and were deployed according to a stratified sampling design.The Alabama coastline was divided into three major sampling areas, each of which contained smaller subareas (Fig. 1).Each year, 240 nets were set for sampling, with specific locations within subareas selected at random with boating conditions allowing.Sampling does not necessarily occur at the exact same coordinates as previous occurrences when a subarea is resampled.Sampling efforts were increased in the summertime, when diversity of catches and variance increases, according to a Neyman allocation; however monthly efforts consistently range from 8 to 13 net sets.Both large and small mesh gillnets were set near each other such that they never intersect or connect, but can be observed simultaneously by sampling crew.Locations of sample sites are shown in Supplementary Fig. S1.A two-tailed t-Test assuming unequal variance indicated that there were no significant differences between annual catch per unit effort (CPUE), calculated as captured individuals/gillnet hours, of bull sharks caught in small, perpendicularset nets vs. large, parallel-set nets (p = 0.31, t = 1.03).Thus, CPUE from both nets was combined for analysis.Sharks were measured to fork length.Environmental conditions including date, time, soak time, tide (high or low), depth, set type, and GPS coordinates were recorded during the soak.Water quality data including salinity (ppt), dissolved oxygen (mg/L), and temperature (°C) were also collected on the surface at the midpoint of the gillnet using a YSI 85.

Analysis
Coastal Change Analysis Program (C-CAP) data was used to obtain LU/LC area coverage for the Mobile Bay watershed, with values available for the years 2001, 2006, 2011, and 2016 38 .These LU/LC values were linked to catch records based on proximate sampling years as temporal resolution of LU/LC allowed (e.g., a shark captured in 2003 would be linked to 2001 LU/LC data, while a shark captured in 2004, would be linked to 2006 LU/LC data).Riverine discharge values were retrieved from USGS gages located at Claiborne and Coffeeville, AL, which mark the upper bounds of the estuary 36 .Riverine discharge values were linked to catch data records based on mean monthly values.Monthly Chl-a values were extracted as raster data (4 km spatial resolution) from the MODIS-Aqua satellite, and linked to catch data by spatial overlap and monthly values.
Time series analysis of annual CPUE were conducted using the 'lm' function, which conducts linear regressions, in the base package 'Stats' in R 39 .Trends were considered significant at p < 0.05.Individual univariate regressions were conducted for the effects of year on mean annual values of in situ (i.e., depth, SST, dissolved oxygen, and salinity) and retrieved (i.e., LU/LC, riverine discharge, and Chl-a) data to identify potential changes in environmental conditions throughout the study period (Table 2).Residuals of the data were assessed for normality using a Shapiro-Wilks test.
Boosted Regression Tree models were fitted for bull shark occurrence (i.e., presence-absence) to develop Species Distribution Models (SDMs) using a total of seventeen variables (Table 1).Boosted regression trees are an additive regression tree model in which models are composed of simple regression trees and built upon an increasing number of iterative trees until the collective model error is minimized and an optimal tree is created.These models are powerful because they are robust to collinearity, missing data, outliers, and can accommodate nonlinear trends 40 .The response variable was binary: presence (1) or absence (0) of a bull shark, therefore the BRTs were built using a Bernoulli distribution and can be interpreted as a probability of presence.Because these SDMs are built on the assumption that environmental conditions at the time and location of a capture indicate suitable habitat, they can be interpreted as habitat suitability models.Boosted regression trees were developed iteratively with the 'gbm.step'function in the 'dismo' package 41 based on different combinations of learning rate (LR), bag fraction (BF), and tree complexity (TC) model parameters 42 .Learning rate describes the proportion of influence for each individual tree term, bag fraction is the proportion of testing data used to train the model, and tree complexity is the number of nodes on an individual tree.Models were fitted in iterations of BF = (0.5, 0.7, 0.9),  input variables in the optimal boosted regression tree model using the 'predict' function in the aforementioned 'dismo' package.This yielded two respective habitat suitability (or probability of occurrence) maps that reflect early and recent portions of the time series.These methods were guided by best practices according to Hijmans and Elith 42 and all analyses were conducted using R version 4.3.1 39 .Using the 'gbm.loop'function in the 'gbm.auto' package 40 , coefficients of variation, representative of model error, were calculated for each raster cell in the habitat suitability maps.

Results
The gillnet dataset included 440 bull sharks caught from 2003 to 2020.Bull shark fork length ranged from 335 to 1068 mm, indicating that all captured individuals were immature 44 .Catch per unit effort significantly increased over five-fold during the 20-year study period from 0.0012 individuals/gillnet hrs to 0.0068 individuals/gillnet hrs (Fig. 2).Concomitantly, mean SST rose from 22.3 °C in 2001 to 23.0 °C in 2020.Overall annual SST increased at a mean rate of 0.05 °C/yr (Table 2), which aligns with the rates observed by Turner et al. 29 .There were also substantial decreases over time in forested (FOR), agriculture (AGR), bare land (BAR), woody wetland (WDW), and Scrub/Shrub (SCB) land cover, and increases in high-intensity development (HID), open-space development (OSD), and low-intensity development (LID) land cover, as well as Chl-a and SST (Table 2, Fig. 3).The final BRT model for bull sharks yielded a CV-AUC of 0.858, suggesting "very good" model fit according to criteria established in Lane et al. 45 .Sea surface temperature accounted for the most variability in suitable bull shark habitat at 27.6%, followed by salinity at 18%.Riverine discharge, Chl-a, and depth were comparable, at 11.6%, 10.6% and 9.4% respectively.SST conveyed the most apparent relationship with model response, with SSTs above ca.22.5 °C resulting in a positive likelihood of bull shark presence (Fig. 4).This likelihood increased in magnitude with warmer SST and exhibited no maximum thermal threshold (Fig. 4).Observable trends are less clear for the other variables.When incorporated into the final model, the random variable yielded an explanatory power of 8.9%.As such, the marginal effects of riverine discharge, Chl-a, and depth indicate they contributed to model performance only somewhat more than random.
Species distribution model results from 2003 indicated greatest habitat suitability in the southwest and southeast portions of the study area as far north as Daphne, AL (Fig. 5a).However, habitat suitability exhibited limited variability across the study area in 2003.By contrast, SDM results from 2020 indicated a wider range of habitat suitability, with least suitable habitat located in the southeast portion of the study area, and most suitable habitat concentrated near Daphne and along the western shoreline of the bay (Fig. 5b).Changes in net predicted relative abundance from 2003-2020 indicated increases in mean habitat suitability that align with the early locations of greatest suitability that expanded in area and magnitude in 2020 (Fig. 6).Mean habitat suitability increased The observed trend in CPUE aligns with the results of the species distribution models yielded by the BRT analysis, which indicate both an increase in the magnitude and geographic extent of habitat suitability for juvenile bull sharks on the Alabama coast over the last 20 years (Fig. 5).Coefficient of variation for model predictions were low (< 1%), and did not appear to exhibit any spatiotemporal patterns (Supplementary Fig. S6a, b).

Discussion
Using long-term gillnet sampling, remote sensing and spatiotemporal modeling, this study documents changes in habitat suitability of juvenile bull sharks in a suspected Gulf of Mexico nursery habitat amidst the changing seascape of the Anthropocene 2 .Results indicate that increasingly warm SST off the coast of Alabama has  coincided with increases in bull shark abundance and habitat suitability, with consequences for future habitat suitability in this region.Given the magnitude of the effects of SST on the SDM, as well as the apparent trend of SST over the study period, we deduced that of the variables tested, increases in SST are primarily responsible for the observed change in juvenile bull shark habitat suitability.This study illuminates the role of increasing temperatures in expanding suitable habitat for bull sharks in the northern Gulf of Mexico.Adaptation occurs at evolutionary time scales, which many forecasters predict will not keep pace with rising temperatures associated with climate change for long-lived species like sharks 46 .However, shifts in animal behavior and distributions can be much more responsive, and may enable mobile species to adjust to changing conditions.Range expansions of the juvenile bull shark have been documented in the northwestern Atlantic as a response to increasing SST, demonstrating the ability of this species to alter their habitat use as such 14 .This study showcases how the use of an estuarine habitat in the northern Gulf of Mexico, a location central to their latitudinal range, has not only persisted, but perhaps been enhanced, in conjunction with warming waters.The lack of an upper SST threshold up to 34.7 °C suggests thermal tolerance for juvenile bull sharks in this area exceeds the thermal maximum of 30 °C postulated by Drymon et al. 23 for the region, but does fall within the range of preferred SST of 27-37.3°C identified for juveniles near southern Florida 47 .These findings also align with those of Curtis et al., who identified no thermal maximum and a preference for water temperatures > 20 °C for juveniles on the Atlantic coast of Florida 48 .Other studies in the Gulf of Mexico found either no significant effect www.nature.com/scientificreports/or a positive with temperatures of 24-26 °C for adult bull shark distributions, both of which contrast with the observed trends for juvenile bull sharks 49,50 .The difference in preferred temperature ranges between juveniles and adults may indicate different physiological constraints in the face of thermal pressures, suggesting the effects of climate change on species may vary based on life history stage.Further, changes in habitat use may also vary ontogenetically, as juveniles may be forced to occupy less suitable habitat to avoid larger conspecifics.Thermal metabolic tolerance of juvenile bull sharks has been found to exceed 30 °C51 , suggesting juvenile thermal tolerance can be robust to warmer temperatures, allowing them to occupy unique habitat.As such, future research into the effects of increasing temperatures on different life and reproductive phases is merited.Contrary to the grim global outlook forecasted for many shark species 2,8 , this study illuminates the resiliency and potential benefits rendered for juvenile bull sharks in the northern Gulf of Mexico in the face of climate change and coastal urbanization.As of the publication date, the bull shark stock in the Gulf of Mexico is unassessed and undergoing its first stock assessment through the Southeast Data, Assessment and Review (SEDAR) process.Our findings are indicative of a local bull shark population that is persisting and perhaps growing, similar to that reported in the western Gulf of Mexico by Froeschke et al. 52 .However, while species distribution studies rely heavily on the assumption that conditions at the time and location of capture reflect suitable habitat for a species, they generally do not consider the condition of the animal.Warming waters have been shown to impose stresses such as increased feeding demands and aerobic metabolism on sharks 53 ; therefore, juvenile bull sharks in Mobile Bay will likely face increasing metabolic demands as waters continue to warm.Increased feeding propensity can be met with increased risk 17 .However, it currently appears that the benefits of increased temperatures, such as decreased development time and expanded habitat suitability, outweigh the costs.In fact, the region in the current study experiencing the greatest increases in habitat suitability also has low hydrological residence times and optimal salinity ranges 54 , conditions which contribute to a highly productive ecosystem 55 .Should a temperature threshold be exceeded though, understanding how populations may respond is of critical importance for management.At present, migrations of bull sharks tagged off Alabama have not surpassed the Florida Keys 50 , and a northward expansion to alleviate thermal pressures would require an extensive expansion of this population's migratory range.To avoid the energetic demands and predation risks associated with extensive migrations, it is possible that bull sharks may spend increased time in riverine habitat in search of cooler waters, conditions which would not necessarily be physiologically limiting for the bull shark over an extended period of time 56 .Continued research into the physiological limits and migratory capacity of this population are warranted.
The recovery and management of elasmobranchs requires a multifaceted approach, which includes enacting spatial protections of critical habitat 57 .Species distribution modeling is a tool that can be used to delineate critical habitat, and has been used to assess the effectiveness of spatial protections for highly migratory species such as bull sharks 18 .Because their predictive power relies upon a set of environmental criteria unique to a species, they can also be used to anticipate changes in distribution in response to hypothetical changes in environmental conditions, such as those induced by climate change.Under the assumption that species relationships to environmental conditions are held constant, they can be used to develop proactive rather than reactive management strategies that have historically characterized shark management 58 .Climate is listed as a significant threat for sharks, and when coupled with pressures from overfishing, could drive many shark species to more extreme categories of the IUCN Red List by 2050 8,59 .In the Gulf of Mexico, the overfishing 60 and invasive species 61 currently threatening coastal ecosystems may be sufficiently exasperated by climate change to impose additional stresses on this subpopulation's habitat suitability prior to the breach of a physiological thermal maximum.As an example of successful sustainable shark fisheries management 11,62 , the United States will need to employ a dynamic and flexible strategy to continue the legacy of success in shark management 6 .Accordingly, these kinds of spatial analyses may be a useful tool to inform fisheries monitoring and management in their task of keeping pace with the pressures of climate change.
Despite the development that has occurred on the Alabama coast over the past several decades 31 , land-use/ land-cover changes were not identified as significant predictors of habitat suitability.Given the increase in juvenile bull shark abundance in the area, it appears that they are not deterred from occupying habitat near urbanized coastline in the northern Gulf of Mexico.This aligns with the findings of Hammerschlag et al. 63 , in which sharks did not actively avoid urban areas near Miami, FL.With increasing abundances of bull sharks near the Alabama coastline, increases in human-wildlife interactions are a likely consequence.Unlike the 'Jaws Effect' that has historically plagued shark management initiatives, surveys of recreational fishermen find that the impacts they perceive sharks will have on fishing opportunities are what primarily hinder their willingness to support shark protection and conservation initiatives 64 .One example of such an impact is through depredation, defined as the "removal of a hooked fish from an angler's line" 64 .Given that bull sharks have been identified as common depredators in the north-central Gulf of Mexico 65 , these changes could hinder stakeholder willingness to support shark conservation and management plans.With a dwindling commercial and increasing recreational shark fishery in the United States, recreational fishermen have the capacity to influence the success and monitoring of shark management plans in the United States through means such as monitoring, engagement in research, enforcement, and advocacy 66 .As such, addressing the concerns of these fishermen and educating them on the significant ecological role sharks play in the health of estuarine systems should be a priority.
To develop a species distribution model (SDM) that included parameters beyond those collected in situ, this study outsourced additional environmental parameters (i.e., Chl-a, riverine discharge, and LU/LC) to develop a www.nature.com/scientificreports/year-round, increased frequency of sampling in warmer months may skew the input rasters and thus the habitat suitability models to more closely reflect the well-mixed, low salinity conditions of Mobile Bay in the summer 36 .Future studies must consider the balance between the benefits of incorporating data sourced from in situ recordings, which offer a more precise understanding of the relationship between species occurrence and environmental parameters when training models, with remotely sensed data, which may reduce bias and provide a more accurate view of overall mean annual conditions when projecting species distributions.Bull sharks depend on nursery habitats for maturation, as they provide ample prey and reduced mortality during an energetically costly and vulnerable life stage 24 .As such, changes in habitat use of juvenile bull sharks in a suspected nursery area like Mobile Bay are likely to have rippling effects on the persistence of a population 24 .The findings of this study suggest that at SST above 22.5 °C, juvenile bull sharks are resilient and may yield a net benefit from large-scale changes wrought by climate change.However, unlike populations in the Atlantic 14 , the Alabama population of bull sharks faces geographic constraints that may cause them to be uniquely vulnerable should a thermal threshold be exceeded.Further, changes in apex predator abundance may have significant ecological and social consequences, altering the structure of the community through predation and risk mediation behavior of prey, while leading to increased human-wildlife interactions 64,68 .Continued monitoring of the local population is imperative for assessing the scope of these potential impacts and when coupled with species distribution modeling, can be used to anticipate changes in distribution due to the impacts of climate change, equipping managers with flexible and proactive means necessary to mediate the impacts of climate change.
https://doi.org/10.1038/s41598-024-54573-0LR = (0.01, 0.005, 0.001, 0. 000001), and TC =(2, 7, 8).The final model was selected to maximize Cross-Validation Area Under the Receiver Operating Curve (CV-AUC), minimize standard error, and maximize Training Data AUC (TD-AUC)44 .R code was guided by Elith & Leathwick43 .After selection of final model parameters, the model was run again to include a random predictor variable, which was populated using the "RAND" function in Excel.The explanatory power of the random variable served as a threshold for predictive performance, and variables yielding values below the predictive power of the random variable were not reported and considered insignificant predictors.Data collected in situ in 2003-2005 and 2018-2020 were rasterized via ordinary Kriging at a resolution of 500 m, which is the approximate width of the narrowest sampling area (Supplementary Figs.S2-S5).This

Figure 1 .
Figure 1.Extent of gillnet survey area along the Alabama coast.Colored areas indicate major sampling areas.Major sampling areas are denoted by number (1-3), and subareas by major area number, then A-E.Inset map denotes extent of sampling coastline in the northern Gulf of Mexico.

Figure 2 .
Figure 2. Catch-per-unit-effort of juvenile bull sharks per year.The dark grey region indicates standard error.

Figure 3 .
Figure 3. Scatterplot of annual environmental parameters with fitted linear regressions.Variables considered to be significant predictors according to BRT output are designated with red regression lines and those considered insignificant with blue regression lines.LU/LC data only have four available values during the study timeframe and represent exact measurements.All other data are mean annual values.The C-CAP categories for LU/ LC include forested (FOR), agriculture (AGR), bare land (BAR), high intensity developed (HID), open space developed (OSD), low intensity developed (LID), woody wetland (WDW), scrub (SCB), open water (WTR), emergent wetland (EMW), and grassland herbaceous (GRS).RD represents riverine discharge.

Figure 4 .
Figure 4. Line plots displaying the marginal effects of predictors performing better than random for bull sharks derived from boosted regression trees.The y-axis is a logit scale, where values of 0 indicate a random likelihood of bull shark presence, values > 0 indicate increased probability, and values < 0 indicate reduced probability.

Figure 5 .
Figure 5. Species distribution model for bull shark derived from boosted regression trees, which predicts suitability of habitat based on environmental data sourced from 2003 (a) and 2020 (b).

Table 1 .
Data source, mean values and ranges for potential environmental predictors used in training boosted regression trees.AMRD is the Alabama Department of Conservation and Natural Resources-Marine Resource Division, which collected data in situ.

Table 2 .
Time series regression results of environmental parameters.Data reflect mean annual changes from the year 2000 to 2020, with the exception of Chl-a sourced from the MODIS satellite, which begins in 2003, and C-CAP data (indicated by *), which begins in 2001 and ends in 2016 and is recorded every 5 years.Variables with missing R 2 values were calculated as generalized least squares models to account for temporal autocorrelation.All other variable metrics are reported from ordinary least squares models.